s01_joined_PCA_plots_971_not_in_LD

started: AL26Apr2019
last updated: AL14Jul2019

Summary

Input sequencing data: 1,850 vars x 715 samples ( 519BC = 258UBC + 260CBC and 197NFE)
Input eigenvectors: 3,219 samples = 715 wecare-nfe + 2,504 kgen

Start_section

Sys.time()

## [1] "2019-07-14 21:08:19 BST"

rm(list=ls())
graphics.off()

library(knitr)

## Warning: package 'knitr' was built under R version 3.5.2

library(dplyr)

## Warning: package 'dplyr' was built under R version 3.5.2

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.5.2

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

base_folder="/Users/alexey/Documents/wecare/ampliseq/v04_ampliseq_nfe/s14_joined_ampliseq_1kg_PCs/s03_all_variants_not_in_LD_971"

opts_knit$set(root.dir = base_folder)

options(stringsAsFactors = F)
options(warnPartialMatchArgs = T, 
        warnPartialMatchAttr = T, 
        warnPartialMatchDollar = T)

#options(error = browser()) # Type Q or c to exit, drop browser level
# https://support.rstudio.com/hc/en-us/articles/200713843?version=1.1.456&mode=desktop
# https://stackoverflow.com/questions/13052522/how-to-leave-the-r-browser-mode-in-the-console-window/13052588 

Read_data

# Sequencing data (for wecare phenotypes: case/control status)
source_folder="/Users/alexey/Documents/wecare/ampliseq/v04_ampliseq_nfe/s12_check_BRCA1_BRCA2_PALB2"
load(paste(source_folder, "s03_exclude_BRCA1_BCRA2_PALB2_carriers.RData", sep="/"))
base_folder="/Users/alexey/Documents/wecare/ampliseq/v04_ampliseq_nfe/s14_joined_ampliseq_1kg_PCs/s03_all_variants_not_in_LD_971"

# Eigenvectors
eigenvectors_file <- paste(base_folder, "ampliseq_1kg_971_3219_100PCs.eigenvec", sep="/")
eigenvectors.df <- read.table(eigenvectors_file, header=T, sep="\t",quote="")

# Eigenvalues
eigenvalues_file <- paste(base_folder, "ampliseq_1kg_971_3219_100PCs.eigenval", sep="/")
eigenvalues.df <- read.table(eigenvalues_file, header=F, sep="\t",quote="")

# 1kg phenotypes (ethnicity and gender)
kg_phenotypes_file <- paste(base_folder, "integrated_call_samples_v3.20130502.ALL.panel", sep="/")
kg_phenotypes.df <- read.table(kg_phenotypes_file, header=T)

# Clean-up
rm(source_folder, eigenvectors_file, eigenvalues_file, kg_phenotypes_file, genotypes.mx, variants.df)

Check data

# List of objects
ls()

## [1] "base_folder"      "eigenvalues.df"   "eigenvectors.df"  "kg_phenotypes.df" "phenotypes.df"

# Expected number of samples in eigenvectors
2504+715

## [1] 3219

# Dimentions of objects
dim(eigenvectors.df)

## [1] 3219  102

dim(eigenvalues.df)

## [1] 100   1

dim(kg_phenotypes.df)

## [1] 2504    4

dim(phenotypes.df)

## [1] 715  24

# Update eigenvectors
rownames(eigenvectors.df) <- eigenvectors.df$FID
eigenvectors.df <- eigenvectors.df[,-1]

Plot eigenvalues

plot(eigenvalues.df$V1, type="b", ylab="Variance",
     main="Top 100 eigenvectors")

plot(eigenvalues.df$V1[1:10], type="b", ylab="Variance",
     main="Top 10 eigenvectors")

rm(eigenvalues.df)

Prepare data for Ampliseq-Kgen plot

Overlapping NFE samples should be removed from Amplise data before plotting

# Make a table with IDs of overlapping NFE
eigenvectors.df[c(3022,3023),c("IID","PC1")]

##                     IID         PC1
## 9_S346_L008 9_S346_L008 -0.01348900
## 2:HG00097     2:HG00097 -0.00844256

nfe_pca <- eigenvectors.df$IID[3023:3219] # re-processed NFE added to Ampliseq
nfe_ampliseq <- sub("2:","",nfe_pca)
nfe.df <- data.frame(nfe_ampliseq, nfe_pca)

# Remove overlapping NFE from ampliseq-kgen eigenvectors
selected_samples <- ! eigenvectors.df$IID %in% nfe.df$nfe_pca
sum(selected_samples)

## [1] 3022

518+2504

## [1] 3022

eigenvectors_ampliseq_kgen.df <- eigenvectors.df[selected_samples,1:6]
"sample" -> colnames(eigenvectors_ampliseq_kgen.df)[1]

# Prepare ampliseq phenotypes 
# (for clarity of PCA plot wecare samples and controls were not separated)
phenotypes.df[c(518:519),c(1,2)]

##                long_ids illumina_id
## 9_S346_L008 9_S346_L008        S346
## HG00097         HG00097        <NA>

phenotypes_ampliseq.df <- phenotypes.df[1:518,c("long_ids","cc")]
table(phenotypes_ampliseq.df$cc)

## 
##   0   1 
## 258 260

"WECARE" -> phenotypes_ampliseq.df$cc[phenotypes_ampliseq.df$cc==1] # This could be named as WECARE-samples
"WECARE" -> phenotypes_ampliseq.df$cc[phenotypes_ampliseq.df$cc==0] # This could be named as WECARE-controls
table(phenotypes_ampliseq.df$cc)

## 
## WECARE 
##    518

c("sample","group") -> colnames(phenotypes_ampliseq.df)

# Prepare kgen phenotypes
phenotypes_kgen.df <- kg_phenotypes.df[,c("sample","super_pop")]
c("sample","group") -> colnames(phenotypes_kgen.df)

# Merge kgen and Ampliseq phenotypes (latest will be on top in the plot)
phenotypes_ampliseq_kgen.df <- rbind(phenotypes_kgen.df,phenotypes_ampliseq.df)
table(phenotypes_ampliseq_kgen.df$group)

## 
##    AFR    AMR    EAS    EUR    SAS WECARE 
##    661    347    504    503    489    518

# Add eigenvectors to phenotypes
eigenphen_ampliseq_kgen.df <- full_join(
  phenotypes_ampliseq_kgen.df, eigenvectors_ampliseq_kgen.df, by="sample")
dim(eigenphen_ampliseq_kgen.df)

## [1] 3022    7

head(eigenphen_ampliseq_kgen.df)

##    sample group         PC1         PC2         PC3          PC4         PC5
## 1 HG00096   EUR -0.01258860 -0.01642160 0.000925435 -0.001809850  0.00780321
## 2 HG00097   EUR -0.00802238 -0.00726931 0.003157220 -0.002825920  0.00838845
## 3 HG00099   EUR -0.01012120 -0.01259730 0.002215210 -0.006756090  0.03568230
## 4 HG00100   EUR -0.00849754 -0.01684940 0.002751330 -0.001556380  0.01011790
## 5 HG00101   EUR -0.01243290 -0.00932382 0.000159626 -0.005677620  0.00907268
## 6 HG00102   EUR -0.00827301 -0.00616160 0.001818280  0.000230261 -0.00955764

tail(eigenphen_ampliseq_kgen.df)

##            sample  group        PC1          PC2          PC3         PC4          PC5
## 3017 95_S517_L008 WECARE -0.0125655  0.000696678  0.000508483  0.00301086  0.000758076
## 3018 96_S236_L007 WECARE -0.0103028 -0.009400760  0.000408240 -0.00126134  0.012597400
## 3019 97_S509_L008 WECARE -0.0116830 -0.016781000  0.003277310 -0.00472534  0.025380700
## 3020 98_S335_L008 WECARE -0.0108617 -0.011973100 -0.000350267 -0.00303147  0.006616230
## 3021 99_S418_L008 WECARE -0.0124556 -0.014528400 -0.002571390 -0.00271848  0.005434640
## 3022  9_S346_L008 WECARE -0.0134890 -0.012876200  0.004023890  0.00388655 -0.003655670

# Clean-up
rm(nfe_pca, nfe_ampliseq, selected_samples, eigenvectors_ampliseq_kgen.df, phenotypes_ampliseq.df,
   phenotypes_ampliseq_kgen.df, nfe.df, phenotypes_kgen.df)

Make Ampliseq-Kgen PCA plot

# Prepare vector fr colour scale
myColours <- c("EUR"="BLUE", "AFR"="YELLOW", "AMR"="GREEN",
               "SAS"="GREY", "EAS"="PINK", 
               "WECARE"="RED")

myColourScale <- scale_colour_manual(values=myColours)

# Static plot
ggplot(eigenphen_ampliseq_kgen.df, aes(PC1,PC2)) + 
  geom_point(aes(col=group)) +
  labs(title="971 variants not in LD", x="PC1", y="PC2") +
  myColourScale

# Interactive plot
plotly_group <- factor(eigenphen_ampliseq_kgen.df$group,
  levels=c("AFR","AMR","EAS","SAS","EUR","WECARE"))

g <- ggplot(eigenphen_ampliseq_kgen.df, aes(PC1,PC2)) + 
  geom_point(aes(col=plotly_group, text=sample)) +
  labs(title="971 variants not in LD", x="PC1", y="PC2") +
  theme(legend.title=element_blank()) + # To suppress the legend title
  myColourScale                         # otherwise it would be "plotly_group"

## Warning: Ignoring unknown aesthetics: text

ggplotly(g, tooltip="text") # By default the tooltip would also show coordinates 

## Warning in dev_fun(file = tempfile(), width = width %||% 640, height = height %||% : partial argument match of 'file' to 'filename'

# Clean-up
rm(myColours, myColourScale, g)

Prepare data for NFE-Kgen plot

The NFE subset in Ampliseq-nfe data had been re-aligned and re-called, so it would be interesting to see
how it is positioned relatively to the same samples aligned/called as in the “original” kgen VCF

# Select samples
eigenvectors.df[c(2504,2505),c("IID","PC1")]

##                     IID         PC1
## NA21144         NA21144 -0.00747477
## 100_S8_L007 100_S8_L007 -0.01326940

eigenvectors.df[c(3022,3023),c("IID","PC1")]

##                     IID         PC1
## 9_S346_L008 9_S346_L008 -0.01348900
## 2:HG00097     2:HG00097 -0.00844256

eigenvectors_nfe_kgen.df <- eigenvectors.df[c(1:2504,3023:3219),1:6]
dim(eigenvectors_nfe_kgen.df)

## [1] 2701    6

2504+197

## [1] 2701

"sample" -> colnames(eigenvectors_nfe_kgen.df)[1]

# Prepare NFE phenotypes 
eigenvectors_nfe_kgen.df[2504:2505,"sample"]

## [1] "NA21144"   "2:HG00097"

sample <- eigenvectors_nfe_kgen.df[2505:2701,"sample"]
group <- rep("Re-processed_NFE",197)
phenotypes_nfe.df <- data.frame(sample,group)

# Prepare kgen phenotypes
phenotypes_kgen.df <- kg_phenotypes.df[,c("sample","super_pop")]
c("sample","group") -> colnames(phenotypes_kgen.df)

# Merge kgen and re-processed NFE phenotypes (latest will be on top in the plot)
phenotypes_nfe_kgen.df <- rbind(phenotypes_kgen.df,phenotypes_nfe.df)
table(phenotypes_nfe_kgen.df$group)

## 
##              AFR              AMR              EAS              EUR Re-processed_NFE              SAS 
##              661              347              504              503              197              489

# Add eigenvectors to phenotypes
eigenphen_nfe_kgen.df <- full_join(
  phenotypes_nfe_kgen.df, eigenvectors_nfe_kgen.df, by="sample")
dim(eigenphen_nfe_kgen.df)

## [1] 2701    7

head(eigenphen_nfe_kgen.df)

##    sample group         PC1         PC2         PC3          PC4         PC5
## 1 HG00096   EUR -0.01258860 -0.01642160 0.000925435 -0.001809850  0.00780321
## 2 HG00097   EUR -0.00802238 -0.00726931 0.003157220 -0.002825920  0.00838845
## 3 HG00099   EUR -0.01012120 -0.01259730 0.002215210 -0.006756090  0.03568230
## 4 HG00100   EUR -0.00849754 -0.01684940 0.002751330 -0.001556380  0.01011790
## 5 HG00101   EUR -0.01243290 -0.00932382 0.000159626 -0.005677620  0.00907268
## 6 HG00102   EUR -0.00827301 -0.00616160 0.001818280  0.000230261 -0.00955764

tail(eigenphen_nfe_kgen.df)

##         sample            group         PC1        PC2          PC3          PC4          PC5
## 2696 2:NA20819 Re-processed_NFE -0.00997535 -0.0161036 -0.002475000  0.006118480 -0.002790070
## 2697 2:NA20821 Re-processed_NFE -0.00942545 -0.0226294  0.690528000 -0.027241300 -0.067820800
## 2698 2:NA20822 Re-processed_NFE -0.01108080 -0.0224294 -0.004567060  0.008997680 -0.010354500
## 2699 2:NA20826 Re-processed_NFE -0.01179990 -0.0229381 -0.008088150  0.011737200 -0.016112000
## 2700 2:NA20828 Re-processed_NFE -0.01296830 -0.0280139 -0.003442760 -0.000977807 -0.000914826
## 2701 2:NA20832 Re-processed_NFE -0.00645055 -0.0133655  0.000427518 -0.001060870 -0.017429200

# Clean-up
rm(sample, group, eigenvectors_nfe_kgen.df, phenotypes_nfe.df,
   phenotypes_nfe_kgen.df, phenotypes_kgen.df)

Make NFE-Kgen PCA plot

# Prepare vector fr colour scale
myColours <- c("EUR"="BLUE", "AFR"="YELLOW", "AMR"="GREEN",
               "SAS"="GREY", "EAS"="PINK", 
               "Re-processed_NFE"="CYAN")

myColourScale <- scale_colour_manual(values=myColours)

# Static plot
ggplot(eigenphen_nfe_kgen.df, aes(PC1,PC2)) + 
  geom_point(aes(col=group)) +
  labs(title="Re-processed NFE", x="PC1", y="PC2") +
  myColourScale

# Interactive plot
plotly_group <- factor(eigenphen_nfe_kgen.df$group,
  levels=c("AFR","AMR","EAS","SAS","EUR","Re-processed_NFE"))

g <- ggplot(eigenphen_nfe_kgen.df, aes(PC1,PC2)) + 
  geom_point(aes(col=plotly_group, text = sample)) +
  labs(title="Re-processed NFE", x="PC1", y="PC2") +
  theme(legend.title=element_blank()) + # To suppress the legend title
  myColourScale                         # otherwise it would be "plotly_group"

## Warning: Ignoring unknown aesthetics: text

ggplotly(g, tooltip="text") # By default the tooltip would also show coordinates 

## Warning in dev_fun(file = tempfile(), width = width %||% 640, height = height %||% : partial argument match of 'file' to 'filename'

# Clean-up
rm(myColours, myColourScale, g)

Final_section

ls()

## [1] "base_folder"                "eigenphen_ampliseq_kgen.df" "eigenphen_nfe_kgen.df"      "eigenvectors.df"            "kg_phenotypes.df"           "phenotypes.df"              "plotly_group"

sessionInfo()

## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.14.5
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] plotly_4.9.0  ggplot2_3.2.0 dplyr_0.8.1   knitr_1.23   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.1        later_0.8.0       pillar_1.4.1      compiler_3.5.1    tools_3.5.1       digest_0.6.19     jsonlite_1.6      evaluate_0.14     tibble_2.1.3      gtable_0.3.0      viridisLite_0.3.0 pkgconfig_2.0.2   rlang_0.3.4       shiny_1.3.2       crosstalk_1.0.0   yaml_2.2.0        xfun_0.7          withr_2.1.2       stringr_1.4.0     httr_1.4.0        htmlwidgets_1.3   grid_3.5.1        tidyselect_0.2.5  glue_1.3.1        data.table_1.12.2 R6_2.4.0          rmarkdown_1.13    purrr_0.3.2       tidyr_0.8.3       magrittr_1.5      promises_1.0.1    scales_1.0.0      htmltools_0.3.6   assertthat_0.2.1  xtable_1.8-4      mime_0.7          colorspace_1.4-1  httpuv_1.5.1      labeling_0.3      stringi_1.4.3     lazyeval_0.2.2    munsell_0.5.0     crayon_1.3.4

Sys.time()

## [1] "2019-07-14 21:08:22 BST"